Высшие симметрии уравнения КдФ

В.Э. Адлер. Классические интегрируемые системы, весенний семестр 2020

Уравнение КдФ

(1)	(Subscript[u, t])=(Subscript[u, 3])+6u (Subscript[u, 1])

допускает бесконечную последовательность высших симметрий нечётного порядка по производным

(2)	(Subscript[u, (Subscript[t, 2k-1])])=(Subscript[f, k])(u,...,(Subscript[u, 2k-1])),   k=1,2...,

то есть, таких уравнений, что [(Subscript[D, t]),(Subscript[D, (Subscript[t, 2k-1])])]=0.

В частности,

(Subscript[u, (Subscript[t, 1])])=(Subscript[f, 1])=(Subscript[u, 1])  – классическая симметрия, отвечающая однопараметрической группе сдвигов x->x+a,
	(Subscript[u, (Subscript[t, 3])])=(Subscript[f, 2])=(Subscript[u, 3])+6u (Subscript[u, 1])  – само уравнение, также классическая симметрия, отвечающая сдвигам t->t+a,
	(Subscript[u, (Subscript[t, 5])])=(Subscript[f, 3])=(Subscript[u, 5])+10u (Subscript[u, 3])+20(Subscript[u, 1])(Subscript[u, 2])+30(u^2)(Subscript[u, 1])  – первая высшая симметрия,
	(Subscript[u, (Subscript[t, 7])])=(Subscript[f, 4])=(Subscript[u, 7])+14u (Subscript[u, 5])+42(Subscript[u, 1])(Subscript[u, 4])+70 (Subscript[u, 2])(Subscript[u, 3])+70(u^2)(Subscript[u, 3])+280u (Subscript[u, 1])(Subscript[u, 2])+70(Subscript[u, 1]^3)+140(u^3)(Subscript[u, 1])  – следующая высшая симметрия,
и т.д..

Линейное пространство конечных линейных комбинаций, порождаемое этой последовательностью как базисом, принято называть «иерархией КдФ». Точнее, это коммутативная и положительная часть иерархии; есть ещё некоммутативная и отрицательная части (ну, и дискретная часть – преобразования Бэклунда), но это мы пока пропустим.

Здесь собрано несколько способов вычисления высших симметрий (всё это будет рассказано за несколько лекций):

◼ исходя из определения. Это трудоёмкий, но зато общий метод, пригодный для любого заданного эволюционного уравнения (то есть, можно либо найти высшую симметрию данного порядка, либо доказать, что её нет), кроме того, он позволяет найти всю алгебру симметрий, включая классические симметрии с зависимостью от x и t;

◼ при помощи оператора рекурсии. Очень эффективный способ, если этот оператор известен. Его вывод будет дан на лекциях, а здесь мы им просто пользуемся;

◼ по явным формулам из квадратичного рекуррентного соотношения. Ещё более эффективный способ, эквивалентный предыдущему;

◼ по явным формулам, с использованием формального обращения преобразования Миуры; это тоже эквивалентно предыдущему;

◼ при помощи мастер-симметрии. У этого метода примерно та же эффективность, что у оператора рекурсии.

Кроме того, имеется метод неопределённых коэффициентов, который вынесен в отдельную программу
06_Homogeneous polynomial equations.nb. Он применим, когда правые части как уравнения (Subscript[u, t])=P, так и симметрии (Subscript[u, T])=Q  –  однородные многочлены. В результате, условие [P,Q]=0 превращается в систему алгебраических уравнений на коэффициенты этих многочленов. Этот метод пригоден также для поиска законов сохранения. Его легко запрограммировать и использовать, в том числе, для поиска новых примеров.

Определим общие функции для символьных вычислений, которые нам потребуются (они уже обсуждались в программе 02_Inversion of the Miura transformation.nb)

(* список динамических переменных в выражении f *)
vars[f_,u_]:=Union[Cases[f,_u,{0,∞}]]

(* df, дифференциал f *)
dif[f_]:=Plus@@(D[f,#]dd[#]&/@vars[f,u])

(* (Subscript[D, x])(f), (Subscript[D, x]^k)(f), полная производная по x от f *)
dx[f_]:=D[f,x]+dif[f]/.dd[u[k_]]:>u[k+1]
dx[f_,k_]:=Nest[dx,f,k]

(* (SubStar[g])(f), оператор линеаризации g, примененный к f *)
star[g_,f_]:=dif[g]/.dd[u[k_]]:>dx[f,k]

(* [f,g]=(SubStar[g])(f)-(SubStar[f])(g) *)
comm[f_,g_]:=Expand[star[g,f]-star[f,g]]

(* E(f), вариационная производная или оператор Эйлера *)
vard[f_]:=Expand[Plus@@((((-1))^(#[[1]]))dx[D[f,#],#[[1]]]&/@vars[f,u])]

(* порядок по производным *)
ord[f_]:=Module[{v=vars[f,u]},If[v=={},-1,v[[-1]][[1]]]]

(* интегрирование по частям *)
int[f_]:=Module[{a=0,b=Expand[f],c,k}, 
While[k=ord[b];     (* вычисляем k - порядок b *)
And[k>0,                   (* выполняем цикл, пока выполняется условие, что k>0 и *)
c=D[b,u[k]];     (* b линейна по старшей переменной *)
Expand[D[c,u[k]]]==0
],
c=Integrate[c,u[k-1]];   (* интегрируем по переменной, следующей по старшинству *)
a=Expand[a+c];                        (* пересчет a и b *)
b=Expand[b-dx[c]]
];
If[ord[b]==-1,    (* если в результате остаток не зависит от u, то это функция от х *)
{Expand[a+Integrate[b,x]],0},
{a,b}
]
]

Тесты

Здесь можно посмотреть, как работают определенные выше функции.

f=u[2]+p[u[1]]

vars[f,u]

dx[f]
dx[f,2]

int[dx[f]]

vard[f]
vard[dx[f]]

◼ Проверка тождества Якоби

f=u[2]+p[u[1]];
g=Sin[u[1]]+u[3];
h=((u[0])^2)+u[3];
comm[f,g]
comm[g,h]
comm[h,f]
comm[f,comm[g,h]]+comm[g,comm[h,f]]+comm[h,comm[f,g]]

◼ Линейные функции всегда коммутируют

comm[u[2]+21u[11]+u[5],u[23]+9u[11]+u[0]]
comm[5u[33]+u[101]+u[500],u[203]+91u[1]-11u[6]]

◼ Уравнения типа Хопфа тоже коммутируют

comm[a[u[0]]u[1],b[u[0]]u[1]]

◼ Симметрия уравнения Бюргерса

ut2=u[2]+2u[0]u[1]
ut3=u[3]+3 u[0]u[2]+3 ((u[1])^2)+3 ((u[0])^2)u[1]
Expand[star[ut2,ut3]]
comm[ut2,ut3]

◼ Подбор коэффициентов в симметрии уравнения Бюргерса методом неопределённых коэффициентов 
(продвинутые обобщения этого примера см. в отдельной программе)

(* уравнение Бюргерса однородно относительно веса w((Subscript[u, j]))=1+j, вес правой части равен 3 *)
ut2=u[2]+2u[0]u[1];

(* составляем многочлен общего вида веса 4, с коэффициентом 1 при (Subscript[u, 3]) *)
ut3=u[3]+a u[0]u[2]+b ((u[1])^2)+c ((u[0])^2)u[1]+d ((u[0])^4);

(* вычисляем коммутатор *)
comm[ut2,ut3]

(* собираем коэффициенты *)
vars[%,u]
CoefficientList[%%,%]
Union[Flatten[%]]

(* приравниваем их 0 и решаем систему *)
Solve[%==0]

(* находим ответ *)
ut3=ut3/.%[[1]]

(* проверка *)
comm[ut2,ut3]

Вычисление симметрии порядка 5 по определению

Задано уравнение (Subscript[u, t])=f. 
Ищем симметрию заданного порядка (Subscript[u, T])=g в общем виде, то есть, g произвольная гладкая функция от t,x, u,(Subscript[u, 1]),...,(Subscript[u, n]). Зависимость от t допускается, при этом

[(Subscript[D, t]),(Subscript[D, T])]=0   ⇔ (∂_t)(g)+(SubStar[g])(f)-(SubStar[f])(g)=0

В этом уравнении следим только за старшей переменной (с максимальным номером). 
Это даёт некоторое несложное уравнение на g.
Решаем его, уточняем вид g. 
Повторяем процедуру, пересчитывая коммутатор, пока g не определится полностью.

Что значит «следим за старшей переменной» – при вычислении коммутатора на бумаге нужно вовремя остановиться и поставить многоточие вместо членов, не содержащих эту переменную. При вычислении на компьютере – не смотрим на всё выражение, а только извлекаем из него какое-то следствие, например, дифференцируя по старшей переменной.

Clear[g]
ut=u[3]+6u[0]u[1];
uT=g[t,x,u[0],u[1],u[2],u[3],u[4],u[5]];
z=Expand[D[uT,t]+comm[ut,uT]];

Если в последней строчке убрать точку с запятой, будет выведено огромное выражение. Нам оно не нужно. Смотрим только, от чего оно зависит и дифференцируем по старшей переменной.

vars[z,u]
D[z,u[7]]

Тут написано -3(Subscript[D, x])((∂_5)(g)) и это должно равняться 0. Значит, (Subscript[u, T])=a(t)(Subscript[u, 5])+g(t,x,u,...,(Subscript[u, 4])). Чтобы экономить буквы, хвост опять обозначаем g, только от меньшего числа переменных. Переопределяем и повторяем все заново.

Clear[a]

uT=a[t]u[5]+g[t,x,u[0],u[1],u[2],u[3],u[4]];

z=Expand[D[uT,t]+comm[ut,uT]];
vars[z,u]
D[z,u[6]]

Имеем (Subscript[D, x])((∂_4)(g))=0, значит, (Subscript[u, T])=a(t)(Subscript[u, 5])+b(t)(Subscript[u, 4])+g(t,x,u,(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])).

Clear[a,b]

uT=a[t]u[5]+b[t]u[4]+g[t,x,u[0],u[1],u[2],u[3]];

z=Expand[D[uT,t]+comm[ut,uT]];
vars[z,u]
Factor[D[z,u[5]]]

Имеем   (Subscript[D, x])((∂_3)(g))=10a (Subscript[u, 1])+a'/3  ⇒ 
(∂_3)(g)=10a u+a'/3x+c(t)   ⇒   
g=10 a u (Subscript[u, 3])+(a'x/3 +c)(Subscript[u, 3])+g(t,x,u,(Subscript[u, 1]),(Subscript[u, 2]))

Clear[a,b,c]

uT=a[t](u[5]+10u[0]u[3])+b[t]u[4]+(a'[t]x/3+c[t])u[3]+g[t,x,u[0],u[1],u[2]];

z=Expand[D[uT,t]+comm[ut,uT]];
vars[z,u]
Factor[D[z,u[4]]]

Имеем   (Subscript[D, x])((∂_2)(g))=8b  (Subscript[u, 1])+20 a (Subscript[u, 2]) +b'/3 ⇒ 
(∂_2)(g)=20a (Subscript[u, 1])+8b u +b'x/3+d(t)   ⇒   
g=20 a (Subscript[u, 1]) (Subscript[u, 2])+8b u (Subscript[u, 2])+(b'x/3+d) (Subscript[u, 2])+g(t,x,u,(Subscript[u, 1]))

Clear[a,b,c,d]

uT=a[t](u[5]+10u[0]u[3]+20 u[1]u[2])+b[t](u[4]+8u[0]u[2])+(a'[t]x/3+c[t])u[3]+(b'[t]x/3+d[t])u[2]+g[t,x,u[0],u[1]];

z=Expand[D[uT,t]+comm[ut,uT]];
vars[z,u]
D[z,u[3]]

Имеем    (Subscript[D, x])((∂_1)(g))=6c (Subscript[u, 1])+60 a u (Subscript[u, 1])+12b  (Subscript[u, 2]) +8a'u/3+2x a'(Subscript[u, 1])+c'/3+x a''/9.

Нетрудно видеть, что это разрешимо, только если a'=0  (на самом деле, если взять a=t, то можно прийти к нелокальной мастер-симметрии из раздела 6, но сейчас мы это не будем рассматривать). Итак, пусть a=const, тогда

(∂_1)(g)=30a (u^2)+12b (Subscript[u, 1]) +6c u +c'x/3+e(t)   ⇒  
 g=30 a (u^2)(Subscript[u, 1]) +6b (Subscript[u, 1]^2)+6c u (Subscript[u, 1])+(c'x/3+e(t)) (Subscript[u, 1])+g(t,x,u)

Clear[a,b,c,d,e]

uT=a(u[5]+10u[0]u[3]+20 u[1]u[2]+30((u[0])^2)u[1])+b[t](u[4]+8u[0]u[2]+6((u[1])^2))+c[t](u[3]+6u[0]u[1])+(b'[t]x/3+d[t])u[2]+(c'[t]x/3+e[t])u[1]+g[t,x,u[0]];

z=Expand[D[uT,t]+comm[ut,uT]];
vars[z,u]
D[z,u[2]]

Имеем  (Subscript[D, x])((∂_0)(g))=4d (Subscript[u, 1])+32 b  u (Subscript[u, 1]) +2b'u +4x b'(Subscript[u, 1])/3+d'/3+x b''/9.

Отсюда получаем b=const и

(∂_0)(g)=4d u+16b (u^2)+d'x/3+f(t)   ⇒ g=2d (u^2) +(16/3)b (u^3)+(d'x/3+f(t))u +g(t,x).

На этом этапе зависимость симметрия от (Subscript[u, n]) определена, остались только коэффициенты. Выведем коммутатор полностью – от него уже почти ничего не осталось, и соберем коэффициенты при u.

Clear[a,b,c,d,e,f]

uT=a(u[5]+10u[0]u[3]+20 u[1]u[2]+30((u[0])^2)u[1])+b(u[4]+8u[0]u[2]+6((u[1])^2)+16/3 ((u[0])^3))+c[t](u[3]+6u[0]u[1])+d[t](u[2]+2((u[0])^2))+(c'[t]x/3+e[t])u[1]+(d'[t]x/3+f[t])u[0]+g[t,x];

z=Collect[D[uT,t]+comm[ut,uT],u[_],Factor]

Отсюда получаем b=d=0,  f=2c'/3,  g=c''x/9+g(t) . Подставим это и пересчитаем.

uT=a(u[5]+10u[0]u[3]+20 u[1]u[2]+30((u[0])^2)u[1])+c[t](u[3]+6u[0]u[1])+(c'[t]x/3+e[t])u[1]+2/3c'[t]u[0]+c''[t]x/9+g[t];

z=Collect[D[uT,t]+comm[ut,uT],u[_],Factor]

Имеем c''=0,  c=(Subscript[c, 1])t+(Subscript[c, 2]),  g=(Subscript[c, 3]), e=6(Subscript[c, 3])t+(Subscript[c, 4]), где (Subscript[c, i]) произвольные постоянные.

uT=a(u[5]+10u[0]u[3]+20 u[1]u[2]+30((u[0])^2)u[1])+(c1 t+c2)(u[3]+6u[0]u[1])+(c1 x/3+6c3 t+с4)u[1]+2/3c1 u[0]+c3;

z=Collect[D[uT,t]+comm[ut,uT],u[_],Factor]

Итак, мы доказали, прямым вычислением, что все симметрии КдФ порядка не выше 5 представляются как линейная комбинация следующих симметрий:

(* высшая симметрия *)
ut5=u[5]+10u[0]u[3]+20 u[1]u[2]+30((u[0])^2)u[1];  

(* само уравнение, симметрия сдвига по t *)
ut3=u[3]+6u[0]u[1];

(* симметрия сдвига по x *)
ut1=u[1];

(* симметрия растяжения *)
utS=3t(u[3]+6u[0]u[1])+x u[1]+2u[0];

(* симметрия преобразования Галилея *)
utG=6t u[1]+1;

(* проверка *)
Expand[D[#,t]+comm[ut,#]]&/@{ut5,ut3,ut1,utS,utG}

Оператор рекурсии

Последовательность (2) можно строить по рекуррентным формулам

(3)	(Subscript[u, (Subscript[t, 2k+1])])=(Subscript[f, k+1])=(R^k)((Subscript[u, 1])),  R=(Subscript[D, x]^2)+4u+2(Subscript[u, 1])(Subscript[D, x]^(-1)),  k=0,1,...

Вычислим симметрии до порядка 25. Для сравнения быстродействия: в программе с неопределенными коэффициентами за 0.5 сек мы досчитывали лишь до симметрии порядка 15.

R[f_]:=Expand[dx[f,2]+4u[0]f+2u[1]int[f][[1]]]

M=13;

Clear[f]
ClearSystemCache[];
Timing[
f[1]=u[1];
Do[f[j]=R[f[j-1]],{j,2,M}]]

Замечание. Функция Timing возвращает список, в котором первый элемент – время вычисления, второй – результат; в нашем случае результат пуст, так как цикл Do ничего не возвращает, только что-то делает внутри себя. Перед обращением к Timing рекомендуется очищать вычисляемую переменную и кэш, так как иначе при повторном запуске ячейки время окажется заниженным – ленивая Mathematica не будет считать заново, а воспользуется сохраненными в памяти результатами (если какая-то мелкая операция выполняется слишком быстро, то для оценки её скорости создается цикл, в котором она повторяется достаточное число раз, причём перед каждым выполнением делается очистка памяти).

Посмотрим на несколько первых членов иерархии (можно заменить 6 на M, и полюбоваться на симметрию порядка 25). Для наглядности, перейдём к индексной записи.

Do[Print["------------------"];
Print[f[j]/.u[n_]:>(Subscript[u, n])],{j,1,6}]

Замечание. В принципе, возможно и сами вычисления делать с использованием индексов, но это не рекомендуется. Лучше делать wrapper – простое преобразование от одной формы к другой, как в приведённом выше примере, и применять её к выводу или вводу. Дело в том, что внутреннее представление выражения u[n] совпадает с ним самим, а для (Subscript[u, n]) оно выглядит несколько неожиданно и это почти всегда приводит к путанице и неприятностям:

FullForm[u[n]]
FullForm[(Subscript[u, n])]

Проверим, что симметрии попарно коммутируют (для больших порядков это занимает много времени, возьмем до M=6)

Table[comm[f[i],f[j]],{i,1,5},{j,i+1,6}]

Квадратичное рекуррентное соотношение

Обозначим (Subscript[f, k])=(Subscript[D, x])((Subscript[a, k])), тогда рекуррентное соотношение (3) переписывается в виде

(4)	(Subscript[a, 1])=u,       (Subscript[D, x])((Subscript[a, k+1]))=(Subscript[D, x]^3)((Subscript[a, k]))+4u(Subscript[D, x])((Subscript[a, k]))+2(Subscript[u, 1])(Subscript[a, k]),     k=1,2,...

Это эквивалентно одному уравнению для производящей функции:

(5)  	(Subscript[a, x x x])+4(λ+u)(Subscript[a, x])+2(Subscript[u, 1])a=0,       a=1/2+(Subscript[a, 1])/(-4λ)+(Subscript[a, 2])/(((-4λ))^2)+...

Его можно один раз проинтегрировать, умножив на 2a, что даёт

(6)	2a (Subscript[a, x x])-(Subscript[a, x]^2)+4(λ+u)(a^2)=λ.

Постоянная интегрирования в правой части определяется по свободному члену (Subscript[a, 0])=1/2. Точнее говоря, к правой части можно добавить произвольный ряд по (λ^(-1)) с постоянными коэффициентами, что отвечает произволу в выборе постоянных интегрирования в (4).  Если сохранять эти постоянные, то просто к (Subscript[a, k]) будет добавляться линейная комбинация (Subscript[a, k-1]),...,(Subscript[a, 0]). 
Полагая постоянные интегрирования равными 0, мы получаем базис из однородных многочленов.

Опять расписывая (6) по степеням λ, мы получаем новое рекуррентное соотношение. В отличие от (3) оно явное (не нужно использовать интегрирование по частям, только дифференцирование и арифметические операции), но зато квадратичное (приходится вычислять длинные суммы по предыдущим членам). В результате скорость вычислений остается примерно такой же. Но, переход от (3) к (6) замечателен тем, что он доказывает то, что все (Subscript[f, k]) являются полными производными и, следовательно, применение (Subscript[D, x]^(-1)) в (3) не встречает препятствий.

Clear[a]
ClearSystemCache[];

Timing[
a[0]=1/2;
Do[a[k+1]=
Expand[
Sum[2a[s]dx[a[k-s],2]-dx[a[s]]dx[a[k-s]],{s,0,k}]-
Sum[a[s]a[k+1-s],{s,1,k}]+
4u[0]Sum[a[s]a[k-s],{s,0,k}]
],{k,0,M-1}]]

(* проверяем, что результаты совпадают *)
Table[Expand[dx[a[k]]-f[k]],{k,1,M}]

Рекуррентные соотношения, основанные на преобразовании Миуры

Уравнение (6) можно переписать в виде следующей системы на ряды v(z), a(λ):

(Subscript[v, x])+(v^2)=(z^2)-u,    a(λ)(v(-z)-v(z))=z,   (z^2)=-λ,

Уравнение для v(z) нам уже знакомо – это преобразовании Миуры (см. программу 02_Inversion of the Miura transformation.nb). Сначала из него находится ряд v(z)=-z+(Subscript[V, 1])/(2z)+(Subscript[V, 2])/(((2z))^2)+..., затем из второго уравнения ряд a. В результате, расчётные формулы не такие громоздкие, как в предыдущем разделе, и считает вдвое быстрее, но, в общем-то, это одно и то же.

Clear[a,V]
ClearSystemCache[];

Timing[
V[1]=u[0];
Do[V[n+1]=Expand[dx[V[n]]+Sum[V[s]V[n-s],{s,1,n-1}]],{n,1,2M-2}];
a[0]=1/2;
Do[a[k]=Expand[V[2k-1]+2Sum[a[k-s]V[2s-1],{s,1,k-1}]],{k,1,M}]
]

(* проверяем, что результаты совпадают *)
Table[Expand[dx[a[k]]-f[k]],{k,1,M}]

Мастер-симметрия

Мастер-симметрией для КдФ называется уравнение (оно получается применением оператора рекурсии к симметрии растяжения, найденной в разделе 1, подробнее это будет на лекциях)

(Subscript[u, τ])=g=x ((Subscript[u, 3])+6u (Subscript[u, 1]))+4(Subscript[u, 2])+8(u^2)+2(Subscript[u, 1])(Subscript[D, x]^(-1))(u).

С его помощью последовательность симметрий (2) строится по рекуррентным формулам

(7)	(Subscript[u, (Subscript[t, 2k+1])])=(Subscript[f, k+1])=(1/(2k-1))[g,(Subscript[f, k])],   k=1,2,...,

где, как обычно,

[g,(Subscript[f, k])]=(Subscript[f, k,*])(g)-(SubStar[g])((Subscript[f, k]))=(∂_τ)((Subscript[f, k]))-(∂_(Subscript[t, 2k-1]))(g).

Обратим внимание, что g содержит, кроме динамических переменных x,u,(Subscript[u, 1]),... еще и нелокальную переменную – интеграл (Subscript[u, -1])=(Subscript[D, x]^(-1))(u). При вычислении коммутатора её нужно дифференцировать по (Subscript[t, 2k-1]). Как это делать, понятно:

(∂_(Subscript[t, 2k-1]))((Subscript[u, -1]))=(∂_(Subscript[t, 2k-1]))(Subscript[D, x]^(-1))(u)=(Subscript[D, x]^(-1))(∂_(Subscript[t, 2k-1]))(u)=(Subscript[D, x]^(-1))((Subscript[f, k]))=(Subscript[a, k]).

Определим расширенную версию функции star для выражений, содержащих (Subscript[u, -1]).

extstar[g_,f_]:=dif[g]/.dd[u[-1]]->int[f][[1]]/.dd[u[k_]]:>dx[f,k]
extcomm[f_,g_]:=Expand[extstar[g,f]-extstar[f,g]]

Задаем правую часть мастер-симметрии и считаем. Да, действительно, получаются высшие симметрии.

g=x (u[3]+6u[0]u[1])+4u[2]+8((u[0])^2)+2u[1]u[-1];

extcomm[g,u[1]]
extcomm[g,%/3]
extcomm[g,%/5]

По скорости этот способ проигрывает оператору рекурсии:

Clear[F]
ClearSystemCache[];

Timing[
F[1]=u[1];
Do[F[k+1]=extcomm[g,F[k]/(2k-1)],{k,1,M-1}]
]

(* проверяем, что результаты совпадают *)
Table[F[k]-f[k],{k,1,M}]